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Langmuir emissions in space are reported to be clumpy and intermittent. The high-frequency wave power 
appears concentrated in spatial packets, be it amidst the solar wind [1] or on auroral field lines [2]. Due 
to the plasma motion relative to the spacecraft, determining the source for the wave free energy in the 
three-dimensional electron distribution function has always been difficult, since the unstable features pass 
by the detector in presumably too short a time to be measured. The range of unstable phase velocities 
and growth rates have generally been estimated rather than determined by unequivocal measurements. 
The analysis of wave-particle interactions in a space environment has taken recently a new turn with 
the development of wave correlators on board rockets [3] and satellites [4]. Such instruments seek to 
identify correlations between the phase of the wave-field and the fluxes of energetic particles. The data 
interpretation is complex, however, and must be backed by a detailed theoretical understanding of the 
wave-particle interaction, including the phase relation for inhomogeneous packets. To this end Langmuir 
packets interact in g with fast electrons can be studied in the appropriate regime by means of particle- 
in-cell simulations, provided that one succeeds in reducing the level of the fluctuations, enhancing the 
signal-to-noise ratio, and incorporating the appropriate boundary conditions. The first results of such 
simulations are presented here as a test and expansion of previous analysis [5]. 

The Langmuir packet dimension L is estimated in the solar wind to be of the order of 10-100 km 
in the direction of the solar wind flow, while on auroral field lines this dimension is of the order of 
100-1000 m along the geomagnetic field. Due to the vastly different plasma densities, in both cases the 
packets contain several tens of wavelengths. When the wave power is concentrated in this manner, the 
resonant function for the traversing electrons is smooth and the homogeneous wave-train singularity at 
the phase velocities is removed [5]. The resulting correlation between the phase of the perturbations in 
the distribution function and the phase of the localized wave determines 

(1) The most unstable phase velocities excited by the streaming electrons, 

(2) The energy exchange between the waves and the electrons, 

(3) The spatial modifications in the downstream non-oscillating distribution function. 

The streaming electrons carry with them information about the wave field. The resistive component 
of the oscillating perturbation is .responsible for the energy exchange with the waves, thereby controlling 
the extent of the packet. When kL < 4w/(fcAu+), where Au+ is the range of velocities with a positive 
slope, the ballistic, linear perturbation contains a resistive component which extends the size of the packet 
downstream. When kL > 4a;/(JkAu + ), the nonlinear component reduces the size of the packet by limiting 
the range of the energy deposition. Thus there exists a natural packet size which is determined by the 
velocity range of the positive slope. Availability of a reliable correlator on board satellites can provide 
an important tool for an analysis of the wave-particle interaction in space, and the numerical as well 
as analytical understanding of the localized Langmuir packet evolution is of great significance for this 
purpose. 

Following the available experimental data in the solar wind [6] and on auroral field lines, our modelled 
electrons consist of a bulk distribution, a power-law tail, and a tenuous beam with a positive gradient 
at large velocities (Figure la). This distribution is used in the theoretical analysis as well as in the 
particle-in-cell simulations. The simulation system resembles the actual physical system by implementing 
non-periodic boundary conditions, which allows injection of realistic distribution functions at the edges of 
the simulation box. The system is thus open, with a prescribed flux of particles that enters from one end, 
interacts self-consistently with the electrostatic field, and leaves at the other end. /This scenario resembles 
the situation of the solar energetic electrons which stream through the solar wind plasma during type III 
bursts and of the auroral electrons which precipitate info the upper ionosphere after having undergone 
modulations in the turbulent region at an altitude of 6000 km. Langmuir waves are characterized by a 
small group velocity; therefore, the resonant electrons cross a packet in a few tens of plasma periods. In 
the simulation we model the wave-packet by a slowly moving, localized, electrostatic pulse v^ith either 
a monochromatic structure or with a small random-phase spread about the unstable wave numbers. It 
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is assumed that the wave was excited locally due to a plasma instability, reaching a weak amplitude of 
e<f)/T£, 0.05, and its finite dimension L parallel to the magnetic field includes tens of wavelengths. We 
ramp-up an external potential for the first 200 steps (covering a time of 24a; -1 ), allowing the electrons 
to adjust to the electrostatic field. When the wave reaches the prescribed amplitude, it interacts self- 
consistently with the plasma. 




Figure 1. (a) Distribution function described by the core, tail, and beam populations 
(extended over five orders of magnitude) and used in the simulations, (b) Dispersion 
relation: u>(k) - solid curve; y(k) - dashed line. The beam destabilizes a narrow range while 
the tail damps a broad range of modes. 

Figure lb shows the results of the linear dispersion relation for the complex frequency (normalized 
to plasma frequency) u = u r + iy vs. wavenumber (normalized to the inverse Debye length) for the 
distribution function of Figure la. The core (c) and the beam (b) distributions are given by stationary 
and the drifting Maxwellians, respectively, 

Fj(v) = (nj/(y/2riij)) exp[— 0.5[(v - Vj)/ Uj ] 2 \ (1) 

where j = 6, c, n&/n c = 5x 10 -4 , 14 = u c > Vb = 10ti c , V e = 0. The tail (t) distribution is described by 

F t {v) = H( v - s)(n t /(V2^u t ))/3[e- a ' 2 (s/v) Q - exp [-0.5(i>/tq) 2 ]] (2) 

where H(x) is the Heaviside function, j3 = e~ a / 2 ^/oc / (2n) / (a - 1) - 0.5[1 - erf(y/aj 2)] is the normal- 
ization factor, u t — u c , n t fn c — 3 x 10 -3 , a is the power-law index, and the junction point with the core 
distribution functions is given by s — ua 1 ! 2 . This choice satisfies the highest possible degree of analyticity 
where the combined distribution F c + F t + Fj, and its first derivative are continuous at v = s. In order 
to evaluate the susceptibility of the tail, we perform a numerical integration over the relevant region 
in velocity space. To ensure the causality of the Vlasov equation, we follow the general procedure of a 
Landau contour [7] by adding the contribution of the deformed contour around the pole for 7 < 0. The 
susceptibilities of the core and of the beam distributions are performed analytically using the standard 
plasma Z functions for the stationary and the drifting Maxwellians, respectively. 

Figure lb shows that the tail electrons damp the waves except for a narrow range of phase velocities 
below t j/k £3 Vi,. The existence of a tail will confine the unstable spectrum to a narrow range unlike 
other simulations of beam-plasma instability [8]. Around the resonant phase velocities a deformation 
of the dispersion relation occurs in both real and imaginary parts due to the beam. This deformation 
disappears only for very dilute beams [9]. For our parameters the instability is kinetic and occurs on the 
Langmuir branch. 

In order to simulate a realistic distribution as given by Eq. 1-2, which covers five orders of magnitude 
in phase-space density, and to resolve the important features at the low-density resonant velocities, 
we divide the macro-particles which occupy the lower density region of phase space into several sub- 
particles with smaller mass and charge values. This procedure, which preserves the charge-to-mass ratio 
and the effective plasma frequency is crucial for meaningful results in our parameter regime. Packet size 
(1 kL ~ 10—100), resolution in k space below the resonant wavenumbers, and a sufficiently large simulation 
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system to analyze the evolution of the packet boundary impose a requirement on a simulation system 
of the order of a few thousand Debye length A u . Sufficient statistics of particles in all regions of phase 
space at the Debye scale-length and reasonably low noise level ( e<f>/T < 0.05) require the total number 
of particles to be greater than one million. 

It has been shown previously [5] that in the linear approximation, for a Gaussian electric field envelope 
E = 2? 0 exp[-(x/2X) 2 ], the perturbed distribution is given by 

f L (x,v,t)= + <* (3) 

where the real argument of the Z function a = kL( l-v p /v) describes the proximity to the phase velocity 
v p of the resonant electrons during their traversal through the packet and the imaginary argument of the 
Z function x = x/2L describes the proximity of the measured distribution from the center of the packet. 
The notable property of the Z function across the real axis gives the asymmetry of the downstream vs. 
upstream distribution. The function fi retains the information of the relative phase and also carries 
its information ballistically to a region outside of the packet. The real part of the complex Z function 
contributes to the reactive part of the perturbed distribution function and affects mainly the group 
velocity of the waves. The imaginary part of the complex Z function contributes to the resistive part 
of the distribution function, which affects the energy exchange between the particles and the wave and 
causes changes in the wave profile. Similarly, the time-independent distribution function which results 
from the deformation of the impinging fluxes along their trajectories is given by 

G(x,v) = ( — ) 7 r^M(v)[e“ 2<r3 (l +er/(x + i<r))(l + erf( X - i<x))] (4) 

where M(v) = kLv p ^-Jfc - v ^ -f and /i = ( 2ekE 0 /mu 2 )(kL ) 2 is the iterative expansion pa- 

rameter. The nonlinear correction which retains the phase information and takes into account the spatial 
deformation of the stationary distribution as well as higher-order coupling of the linear perturbation is 
given by 

= xfEp** k '- ut) '-t-I-„dy e -3yV-( y -*)£ [(^) 4 (P(t/, u) + iQ(y, t,))] + cc (5) 

where P(*,u) = v p ^ [Z r (v) -a \ Z{u ) | 2 ] - - v^) I Z(v) | 2 and 

Q{x, v) = VpT? [Z{v) ~ V2Z(V2v)} + ^ Z 2 {u) for u = -a - i X . 

The resistive part of the linear perturbation and the nonlinear perturbation g constitute the 
source of the energy exchange between the waves and the particles. For >0.6 the contribution of 
the nonlinear perturbation becomes comparable to that of the linear perturbation. Figure 2 shows the 
calculated contribution of the peTturbed current and the energy deposition to the packet due to both of 
the terms. The spatial profile of the linear contribution approximately follows the shape of the packet, 
resulting in a self-similar growth of the field profile. On the other hand, the non-linear contribution has 
an opposite sign and its profile is shifted downstream relative to the packet position. The combined 
contribution is therefore asymmetric with respect to the wave profile, enhancing its amplitude upstream 
and shrinking its size downstream. 

Contributions to the current 
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Figure 2. Contributions of the linear (upper solid line), the nonlinear (lower solid line), and 
the total (dashed line) current, (in arbitrary units) to the energy deposition, as a function 
of the spatial coordinate. The total contribution indicates a localization of the packet. 
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Figure 3 shows several snapshots of the wave profile as recorded from a simulation run. The original 
size of the packet is kL — 25 ; it is ramped-up adiabatically for the time period of 24a> -1 , and afterwards 
the packet evolves self- consistently. For low values of fi < 1, the electrons are perturbed according to Eq 4; 
they exhibit the deformation during the packet traversal and carry the ballistic information downstream. 
As the packet amplitude is growing, the fresh electrons which enter the system are affected by the 
stronger electric field. As a result their stationary deformation G increases, and the phase information in 
the nonlinear perturbed distribution function becomes spatially dependent. The resulting effect on the 
packet is its growth and localization, as shown at a later time of the simulation. 
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Figure 3. Snapshots of the potential $(:c) in units of eT as a function of x (in units of A^>), 
recorded at different times of the simulation. 
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